close all;
clc, clear;
warning off;
%% 读取数据
X = xlsread('data4.xlsx');
X = X(:, 2:end); % 特征数据
Y1 = xlsread('newdata.xlsx');
Y0 = Y1(:, 2); % 目标变量
%% 输入改变
ct=10;%12,10
numi=600;%800,600
study=0.1;%0.01,0.1
numstudy=numi-200;
data=[1038547,1032736,1049427.75000000,1049677,1058578.37500000,1050872.75000000,1054769.25000000,1075204.75000000,1075717.75000000,1080017];
%% 数据归一化
[x,psin]=mapminmax(X',0,1); 
[y,pcout]=mapminmax(Y0',0,1);
num = length(Y0); % 总样本数
state=randperm(num);%打乱样本
r = 0.7; % 训练集占比
trainnum = floor(num * r); % 训练集样本数
xtrain=x(:,state(1:trainnum));
ytrain=y(state(1:trainnum))';
%% 适应格式
for i=1:trainnum
    x_train{i,1}=xtrain(:,i);
end
%% 创建模型
layers=[
    sequenceInputLayer(5,"Name",'input')             %输入层
    lstmLayer(ct,'OutputMode','last')        %LSTM层
    reluLayer                               %Relu激活层
    fullyConnectedLayer(1)                  %全连接层
    regressionLayer ];                      %回归层
%% 参数设置
options=trainingOptions('adam' ,...                              %梯度下降算法
                        'MaxEpochs',numi, ...                     %最大迭代次数
                        'InitialLearnRate',study, ...             %初始学习率
                        'LearnRateSchedule','piecewise', ...     %学习率下降
                        'LearnRateDropFactor',0.1,...            %学习率下降因子
                        'LearnRateDropPeriod',numstudy, ...           %经过多少次训练后学习率为*0.1
                        'Shuffle','every-epoch', ...             %每次训练打乱数据集
                        'Plots','training-progress', ...         %画出曲线
                        'Verbose',true);           
%% 训练LSTM模型
[net,trainInfo]=trainNetwork(x_train,ytrain,layers,options);
%% 仿真预测
re1=predict(net,x_train);
%% 未来预测（10年）
futurenum = 10; % 预测未来10年
outputs = zeros(futurenum, 1); % 初始化预测输出
x=x';
rate = mean(diff(x(end-20:end, :))); % 使用最近20年的数据计算平均增长率
% 初始化未来X矩阵
XFuture = zeros(futurenum, size(x, 2));
% 线性增长的特征
for i = [1, 2, 3, 5]  % 对第1、2、3、5列特征，假设是线性增长的
    Xpt = x(:, i);
    years = (1:size(x, 1))';
    p = polyfit(years, Xpt, 1); % 线性拟合
    fyears = (size(x, 1)+1:size(x, 1)+futurenum)';
    lgrowth = polyval(p, fyears);
    % 添加更大的随机波动
    XFuture(:, i) = lgrowth + randn(futurenum, 1) * 0.05;
end
% 平稳的特征
XFuture(:, 4) = repmat(x(end, 4), futurenum, 1) + randn(futurenum, 1) * 0.01;
for i = 1:futurenum
    XFuture1=XFuture(i,:)';
    % 预测未来的Y值
    output = predict(net, XFuture1);
    outputs(i) = mapminmax('reverse', output, pcout); % 反归一化
end
disp('当前指标的波动率')
p=sum(outputs'./data)/10;
disp(num2str(1-p));

